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NUMERICAL ANALYSIS AND PARAMETRIC STUDIES OF THE BUCKLING 
OF COMPOSITE ORTHOTROPIC COMPRESSION AND SHEAR PANELS 


Jerrold M. Housner and Manuel Stein 
Langley Research Center 

SUMMARY 

A computer program has been developed for the combined compression and shear 
of stiffened variable thickness orthotropic composite panels on discrete springs; boundary 
conditions are general and include elastic boundary restraints. Buckling solutions are 
obtained by using a newly developed trigonometric finite -difference procedure which 
improves the solution convergence rate over conventional finite-difference methods. The 
trigonometric finite-difference procedure introduces two new parameters into the solution. 
These parameters can be computed by the program or selected by the user. The validity 
of the program has been substantiated by eomparisons with existing solutions, and a pro- 
gram listing, input description, and sample problem are provided. 

The classical general shear -buckling results (in terms of universal orthotropic 
parameters), which exist only for simply supported panels over a limited range of ortho- 
tropic properties, have been extended to the complete range of these properties for sim- 
ply supported panels and, in addition, to the complete range of orthotropic prqjerties for 
clamped panels. The program has also been applied to parametric studies which exam- 
ine the effect of filament orientation upon the buckling of graphite -epoxy panels. These 
studies included an examination of the filament orientations which yield maximum shear 
or compressive buckling strength for panels having all four edges simply supported or 
clamped over a wide range of aspect ratios. Panels with such orientations had higher 
buckling loads than comparable, equal -weight, thin-skinned aluminum panels. Also 
included among the parameter studies were examinations of combined axial compression 
and shear buckling and examinations of panels with rotational elastic-edge restraints. 

INTRODUCTION 

The use of filamentary composite materials in aircraft and space structures offers 
a potential for weight savings over conventional (all metal) construction. Also, compos- 
ites introduce added versatility into the design process by allowing the structure to be 
better tailored to meet the design criteria. One such design criterion is the prevention 
of compressive and shear buckling in panels of laminated construction. In laminated 



panels the stiffness properties can be tailored by controlling the filament orientation in 
each lamina, 

A considerable amount of literature exists on the buckling of flat isotropic and 
orthotropic panels under various boundary conditions. (See refs. 1 to 6.) Few results 
exist, however, for finite aspect -ratio panels, especially for shear buckling of orthotropic 
panels. General results for shear buckling, in terms of universal orthotropic parameters, 
exist only for simply supported panels over a limited range of orthotropic parameters. 

(See ref. 6.) Several general-purpose computer programs exist which could be employed 
to obtain results for panels with general boundary conditions under general loading states 
(refs. 7 to 9). These programs, however, tend to be expensive to use in performing 
parameter studies; therefore, a program which is suitable for performing parametric 
buckling studies of orthotropic flat rectangular panels was developed and is employed in 
this paper. 

The present computerized analysis is applicable to the combined compression and 
shear buckling of stiffened, variable -thickness, flat rectangular orthotropic panels on 
discrete springs; boundary conditions are general and include elastic boundary restraints. 
Calculation of the flexural stiffnesses of a laminate from the properties of filament - 
reinforced laminas is automatically performed. The analysis makes use of a newly 
developed trigonometric finite -difference procedure. In contrast to conventional (poly- 
nomial) finite differences, trigonometric differences take advantage of the sinusoidal 
form of the buckle pattern to achieve converged solutions with fewer degrees of freedom, 
hence reducing computer time. The analysis has been validated by many comparisons 
with solutions in the literature and has been used to produce a variety of additional ortho - 
tropic and some isotropic panel results. 

The classical general results for the shear buckling of simply supported orthotropic 
panels are extended in this paper to cover the complete range of orthotropic parameters. 
Also, the general results for the shear buckling of clamped panels over the complete 
rai^e of orthotropic parameters have been calculated and are presented herein. In addi- 
tion, it is of practical interest to present results which consider the effects of filament 
orientation i:^on the buckling strength of laminated composite panels. Consequently, 
parameter studies are presented for graphite -epoxy panels of various aspect ratios, 
boundary conditions, and in -plane loadings over a wide range of f ilament orientations, and 
those orientations which led to maximum buckling loads are identified. Finally, results 
are presented for the shear buckling of simply supported isotropic panels, each with a 
central stiffener. 
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SYMBOLS 


a,b dimensions of rectangular plate parallel to X- and Y-axes, respectively 

coefficients defined by equation (C3) 

Cx,Cyx correction factors defined in equations (B4) and (B5) 

D isotropic plate flexural stiffness 



El flexural stiffness of discrete stiffener 

Ej^,E 2 Young's moduli of fibrous reinforced material parallel to fibers and trans- 
verse to fibers, respectively 

Gj 2 shear modulus of fibrous reinforced material 


h core thickness of sandwich plate 

Il,l 3 row designations of boundaries (D and (3) (see fig. 2(a)) 

J2?^4 column designations of boundaries (D and (3) (see fig. 2(a)) 

k| discrete lateral spring stiffness 

kj^ uniformly distributed rotational spring stiffness 

kg shear ^buckling load coefficient 

kx.% stiffness of rotational springs which resist moments acting about Y- and 

X-axes, respectively 



3 



plate stiffness terms defined by equation (A13) 


M,N total number of rows and columns of finite -difference stations, respectively 

Me,Ng total number of rows and columns of finite -difference stations at which equilib- 

rium is satisfied 

Mx,My,Mxy bending moments in plate (see fig. 1) 


Nx,Ny,Nxy 


in-plane loads (see fig. 1) 


Nx,Ny,Nxy 


shear-buckling stress coefficients 
respectively 


b^Nx 


b^Ny b^Nxy 


^xj^yj^xy 


buckling parameters 
respectively 




b^N- 


xy 


Eit‘ 


i-m 


N- 


X ,N- 

Xq> . 


xyo 


P 

^x>^y>^xy 

I^xj^xy 


t 

txjtyjtxy 


w 


buckling loads for pure axial compression and pure shear, respectively 
buckling eigenvalue (see eq. (13)) 

change of Nx, Ny, Nxy with p, respectively (see eq. (13)) 
ratio of Nx/Nxq and Nxj^Nxy^) respectively 

spring-stiffness terms defined by equation (A17) 

total thickness of sandwich plate 
values of l^x> ^^y^ ^^xy v^hen p * 0 

displacement of panel in positive z -direction 


x,y,z 


panel coordinates shown in figure 1 



curvature terms defined in equation (A14) 



6U internal virtual work 

6Vj^ virtual work of in -plane loads 

6Vg virtual work of discrete springs 

Ax, Ay f inite -difference mesh spacings in x- and y-directions, respectively 

A ^ 

Ax,Ay trigonometric finite -difference coefficients as defined in equation (10) 

Ax, Ay trigonometric finite -difference terms defined in equation (A24) 

9 filament orientation (see fig. 2(a)) 

©,B universal orthotropic parameters defined in equations (15) and (16) 

\x)^y trigonometric parameters defined through equation (10) 

v^2 major Poisson ratio relating contraction normal to filament direction to 

extension parallel to filament direction 

Ixjlyj^x’^^y functions defined by equations (A6) to (A9) 

Xij twist terms defined in equation (A16) 

iji'y curvature terms defined in equation (A15) 

Comma preceding a subscript denotes differentiation with respect to the subscript. 
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ANALYSIS 


Assumptions 

The buckling analysis of linear elastic orthotropic plates has been carried out under 
the following assumptions: 

1. Coupling between bending and extensional deformation is neglected. (In practice 
this assumption implies a midplane symmetric laminated panel.) 

2. Coupling between bending and twisting deformation is neglected, (In practice 
this assumption implies a balanced laminate.) 

3. The deformations of the panel obey the Kirchhoff hypothesis (see ref. 10). 

4. The nonlinear strain -displacement relationships used to obtain (linear) buckling 
equations are 

ex = e X + !('>', x)^ 

ey = ''.y iKy)^ 

’'xy = “,y + V + w,xW y 

where ex, Oy, and yxy the strains and u, v, and w are the displacements in 
X-, y-, and z -directions, respectively. 

5. The in -plane loads, Nx, Ny, and Nxy, are unif ormly distributed along the 
appropriate edges of the plate, 

6. Discrete stiffeners have no torsional stiffness and are symmetrically disposed 
with respect to the neutral surface of the panel. 

Governing Equations 

The internal virtual work of the panel during buckling may be expressed as 


pb pa 

5U = Jq (^x^w^xx + ^y^^,yy + 2Mxy6w^xy) ^ 


( 1 ) 


where a and b are the dimensions of the panel parallel to the X- and Y-axes, respec- 
tively, and 6 is the variational operator. Also, 
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Mx = Diiw^xx + Di2W,yy'^ 

My = Di2W^xx + ^22W^yy > 

Mxy = 2D00W^xy 


( 2 ) 


The sign conventions of the bending moments are given in figure 1, and the flexural stiff- 
nesses, Dji, Dj2, D 22 , and D 00 , given in reference 11, are about a unique neutral 
plane which has the property that matrix [B], which represents coupling between bending 
and extension, is null with respect to this plane. As given by reference 12, the virtual 
work of the applied in -plane loads is given by 

^ ^ ^NxW x^w X NyW y^^,y ^xy^,y^^,x ^xy^,x^^,5^ 

where the sign conventions for Nx, Ny, and Nxy are shown in figure 1. 

In appendix A, equations (Al) to (A3) are expressed in trigonometric finite -difference 
form (see fig. 2 for finite -difference station layout) and are substituted into the statement 
of the principle of virtual work, that is. 


6U = 6Vn + 6Vg (4) 

where 6Vg is the virtual work of the discrete springs. (See appendix A, eq. (All).) 
Equation (4) yields the governing equations which are of the following form: 


Kij + Sij 




%^ij ^^xyXij 0 



(5) 


where Ky, Sy, oty, and Xy are defined by equations (A13) to (A17) in 
appendix A. 

The nvimerical technique of trigonometric finite differences and the numerical 
extraction of the buckling loads Nx, Ny^ and Nxy fi'om equation (5) are different from 
those conventionally used and therefore require further discussion. 


Numerical Techniques 

Trigonometric finite differences . - Conventionally, the central difference approxi- 
mation for the derivative of a function f(x) at x = xq is approximated as 
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( 6 ) 


1^0) 


f Xn + 





The right-hand side of equation (6) is denoted as the conventional finite-difference approx- 
imation for the derivative. In the limit as the finite -difference mesh spacing Ax 
approaches zero, the right-hand side of equation (6) expresses the definition of the deriv- 
ative. If f(x) is parabolic in the neighborhood of xq, 


f(x) = + 72^^ " ^O) + rzi^ - (’7) 

and it may be readily shown that the approximate expression given by equation (6) becomes 
an equality. If, however, f(x) is trigonometric about x = Xq, 


where 


^x 


,, , . 77(X - Xq) 77(X - Xq) 

f(x) = y, + Sin — + y„ COS 

Ax Ax 


is a wavelength parameter. It may be readily shown that 


( 8 ) 


dx 


(xq) = ^ 


■‘x 


f Xn + ■ 



(9) 


where 



7T 


2Ax 



( 10 ) 


The right-hand side of equation (9) is denoted as the trigonometric finite-difference 
approximation for the derivative. ^In a two-dimensional problem a similar set of rela- 
tionships would be derived for the y-direction, introducing the quantities Ay, Ay, 
and Ay.^ 


The only difference between the right-hand side of equation (9) and that of equa- 
tion (6) is that in the trigonometric expression i/Ax replaces 1/Ax of the conventional 
expression. As Ax approaches infinity. Ax approaches Ax and, consequently, the 
trigonometric difference expression reduces to the conventional expression. 
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Convergence of trigonometric finite-difference solutions .- Inasmuch as the buckling 
mode shape is usually trigonometric in nature, the trigonometric finite-difference solu- 
tion can be made to exhibit a much faster convergence rate than the conventional differ- 
ence solution by appropriate selection of and Xy, This advantage is demonstrated 
with several isotropic plate examples discussed in appendix B. The convergence rate 
can also be degraded, however, by an inappropriate choice of X^ and Xy. It should be 
emphasized though, that the selection of ^x and Xy does not constrain the buckle mode 
shape to have wavelengths given by ^x and Xy . Rather , the trigonometric solution will 
always converge to the exact solution if enough degrees of freedom (finite-difference sta- 
tions) are used. 

Selection of trigonometric parameters Xx and Xy.- Selecting appropriate values 

of Xx and Xy which improve the convergence rate of solutions is predominantly based 
upon engineering judgment and experience. One engineering approach which has proven 
useful is to select Xx and Xy based upon the buckle length of infinitely long panels; 
that is. 


^ _ bM 

a /3 


( 11 ) 


Xy 

^ 1 ( 12 ) 

b 

where /3 is the wavelength parameter of an infinitely long panel, defined as the ratio of 
the panel width to the buckle length. The value of /3 for the combined compression and 
shear bucklipg of simply supported and clamped infinite panels may be determined from 
equations (B2) and (B3) in appendix B. Additional suggestions for the selection of Xx 
and Xy are given in appendix B. 

Stability determinant evaluation and eigenvalue extraction .- In this analysis the 
order of the stability determinant is kept to a manageable size by using the two- 
dimensional marching procedure outlined in appendix C. This procedure is basically an 
extension of the one-dimensional procedure used in reference 13. Briefly, the marching 
procedure successively operates on the equilibrium equations at each finite-difference 
station to achieve a relatively low-order stability determinant. 

In searching for the combined load system which produces buckling, it is convenient 
to introduce dimensionless stress coefficients, Nx, Ny, and Nxy, which may be deter- 
mined from the dimensional quantities, Nx, Ny, and Nxy (fig. 1), by multiplying by the 

factor b^ 7 r/Dii. It is assumed that Nx, Ny, and Nxy are linear functions of an eigen- 
value p, that is, 
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Nx = ^ + Pi-x 

Ny = ty + pry > 
Nxy = txy + prxy^ 


( 13 ) 


This assumption allows some loads to be held constant while others are increased to 
buckling, or it allows the loads to increase with a fixed proportionality. 

To find the lowest value of p which makes the stability determinant vanish, a 
determinant plotting technique is used, In order to increase the speed of the plotting 
technique, a variable step size is employed. This step size is based upon a numerical 
parabolic extrapolation of the stability determinant at each step of the determinant plot- 
ting procedure. 


COMPUTEE PROGRAM 


A computer program denoted BOP (Buckling of Orthotropic Panels) has been devel- 
oped for the buckling of flat rectangular orthotropic laminated panels. The program is 
applicable to panels with compression and/or shear loading, discrete lateral deflection 
and rotational springs, discrete stiffeners, and general boundary conditions. 

The program utilizes trigonometric finite differences to improve the problem con- 
vergence and thus requires the selection of and Ay. The user has the option of 
determining and supplying Ax and Ay (based upon the discussion in appendix B) or 
allowing the program to automatically calculate and use values based on equations (11) 
and (12). 

In addition, the user has the option of either (1) Supplying the bending stiffnesses of 
the panel or (2) supplying the elastic moduli, filament orientation, and thickness of each 
lamina in a laminated panel and allowing the program to calculate the bending stiffnesses. 
When the second option is chosen, the program prints the flexural stiffness matrix D, 
defined in reference 11, as well as the laminate Young's moduli, shear modulus, and 
Poisson's ratios. (The second option may be used independently of the buckling analysis.) 
A complete description of the program is provided in appendix D. 

Results from the computer program have been compared with many classical 
results for unstiffened isotropic and orthotropic panels under various boundary conditions 
and with some classical results for stiffened isotropic panels. These comparisons which 
are discussed in subsequent sections were found to be excellent, thereby indicating the 
validity of the program. 
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RESULTS AND DISCUSSION 


Shear Buckling of General Orthotropic Panels 

From the general fourth-order equation for the shear buckling of orthotropic panels 
the buckling load coefficient may be expressed as 


b^N. 


ks - 


xy 


AP11E>22 


(14) 


This coefficient is a function of only two variables 

Q _ ^^ 11^2 2 

DS 


and 

b\4j^ 

- \l5^ 


(15) 


(16) 


where D 3 = D^2 2 D 00 . (Note that an isotropic panel implies 0 = 1.) 

Classically, general shear-buckling results for simply supported finite aspect-ratio 
panels have been obtained only for values of © = 1 (see ref. 6 ). In figure 3 numerical 
results for 0 < 1 have been presented. Also, for completeness and comparison pur- 
poses numerical results for 0=1 are presented. The good agreement between these 
curves and those of reference 6 indicates the validity of the numerical results from the 
computer program. General results for the shear buckling of clamped panels, further- 
more, do not appear in the literature for any range of 0 with the exception of 0 = 1 
(the isotropic case) ; consequently , numerical results for clamped panels are presented 
in figure 4. 

Both the results for simply supported and elamped panels indicate that the percent- 
age decline in buckling load from B = 1 to B = 0 decreases as © increases. Also, 
a comparison of figures 3 and 4 shows that the percentage increase in buckling load of 
clamped panels over simply supported panels increases with increasing 0. The abrupt 
changes in slope appearing in these figures are due to changes in mode shape (from sym- 
metric to antisymmetric modes). As anticipated from isotropic results (ref. 1), these 
abrupt changes are more predominant in clamped panels than in simply supported panels. 

Tables 1 and 2 present the shear -buckling load coefficients used in obtaining the 
general orthotropic panel results of figures 3 and 4. Additionally, the trigonometric dif- 
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ference parameters ^the mesh-spacing parameters a/Ax and b/Ay and the wavelength 
parameters A.x/a and ^y/b) used in obtaining the buckling coefficients are presented 
in tables 1 and 2. 


Shear Buckling of a Simply Supported Panel With a Central Stiffener 

Figure 5 presents results for the shear buckling of simply supported isotropic 
panels each of which contains one central flexural stiffener parallel to either the lor^er 
or shorter edges of the panel. As anticipated, the use of a central stiffener always pro- 
vides an increase in the shear -buckling stress coefficient over that of the imstiffened 

panels = oj . The percentage increase over unstiffened panels is greater in square 

panels than in rectangular panels. In rectangular panels of the same aspect ratio, the 
percentage increase over unstiffened panels is greater when the stiffeners are parallel 
to the longer direction than when they are parallel to the shorter direction. The central- 
stiffener results of figure 5, moreover, are in reasonably good agreement with similar 
results given in reference 14 for slightly curved panels. This agreement indicates the 
validity of the computer program for the solution of stiffened panels. 


Parametric Studies of Orthotropic Filament Reinforced Panels 

Results are presented for the buckling of sandwich panels whose upper and lower 
skins are of laminated graphite- epoxy construction. Although some of the results in this 
section could be obtained from general orthotropic curves, such as those of figures 3 
and 4, it is of Interest to examine the effect of filament orientation upon the buckling load. 
(The material properties for the graphite -epoxy skins are given in table 3, with their 
equivalent general orthotropic parameter values 0 and B at various filament 
orientations.) 

In addition to the assumptions listed in the analysis section of this report, it is 
assumed in this section that 

1. The panel is symmetric about the middle surface 

2. Each lamina has the same filament orientation 6 except for sign 

3. The core carries no load and undergoes no transverse shear deformation 

As a consequence of these assumptions, it may be shown that the buckling param- 
eters Nx, Ny, and Nxy defined as 


Nx = 


b%x 



(17a) 
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(17b) 



(17c) 


depend only on the magnitude of 9 , the panel aspect ratio, and the boundary conditions. 
They do not depend on the thickness of each lamina, the number of laminas, or the core 
thickness. However, in order for assumption 2 of the analysis section to be reasonable — 
that is, neglect of bending-twisting coupling — it may be necessary that the ratio of core 
thickness to total thickness h/t be nearly unity and that the amount of material in either 
cover oriented in the +6 and -6 directions be equal. 

The variation of the buckling load with filament Orientation for panels of various 
aspect ratios is presented in figure 6 for axial compression and in figure 7 for shear. 

The figures indicate that the buckling loads are highly dependent upon filament orientation 
and that optimum orientations (those which yield a maximum buckling load) may be deter- 
mined for each aspect ratio. Also, the figures indicate that clamping has a greater effect 
on compressive buckling than on shear buckling. 

An indication of the buckling strength of the epoxy panels as compared to equal- 
weight aluminum panels is provided by a comparison of the discrete buckling loads 
appearing on the right-hand ordinate of figures 6 and 7 with the curves in the same fig- 
ures. These comparable values are valid for thin-skinned sandwich panels which have 
the same core, of thickness h, as the graphite -epoxy panels, but which have aluminum 
skins. For all the cases considered, a range of filament orientations exists for which 
the buckling strength of the graphite -epoxy panels exceeds that of the comparable alumi- 
num panel with the same aspect ratio and boundary conditions. In the case of a clamped 
square panel in shear, the buckling strength of the graphite -epoxy panel exceeds that of 
the aluminum panel at all filament orientations. 

It should be noted that, if the restriction that each lamina have the same filament 
orientation ±9 is removed, isotropic skins can be produced from groups of three or 
more laminas (for example, 0, +60, and -60) which will have the same weight as the ±9 
skins but will yield a higher buckling load for each case shown in figures 6 and 7 and for 
many other shear and compression loadings. However, this is not necessarily true in all 
cases; for example, in the transverse compression of long panels (a/b approaching zero), 
an orthotropic panel with filaments running transverseiy {9 = 0°) provides a higher 
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buckling load than an equivalent isotropic panel. Furthermore, there are many applica- 
tions where for various reasons (for example, strength or fabrication criteria) ortho - 
tropic panels are preferable to isotropic ones. 

In figures 8 to 11 optimum filament orientations are shown for all aspect ratios. 

The curve of figure 8 was determined from the exact closed-form relationship for the 
compression of simply supported plates (ref. 6), while the curves of figures 9 to 11 were 
determined using program BOP. The abrupt changes in the slopes of these curves are 
caused by changes in the buckling mode shape associated with the optimum filament ori- 
entation. Except for figure 8, the location of these abrupt changes has been approximated 
since it is difficult to determine exactly where they occur. 

In the compressive buckling curves (figs. 8 and 9) the optimum filament orientation 
for small aspect ratio a/b is 0° (parallel to the X-axis or to the direction of compres- 
sion). This orientation angle rapidly increases at about a/b = 0.56 for simply supported 
panels and at about a/b = 1.05 for clamped panels. However, a eomparison of the 
aspect -ratio 1 and 1.1 curves for a clamped panel as shown in figure 6 indicates that the 
optimum buckling load does not exhibit such a rapid charge but decreases slightly as the 
aspect ratio goes from 1 to 1.1. For higher aspect ratios the optimum orientation oscil- 
lates with decreasing excursion about ±45<^ and, in general, a practical filament orienta- 
tion for a/b > 1 is 0 = ±45°. 

In the case of shear buckling (figs. 10 and 11), the symmetry of the problem requires 
that the deviation of the optimum filament orientation from 45° for a panel of aspect ratio 
a/b be equal but opposite to that of a panel with aspect ratio b/a. Also, the peaks of 
figure 7 are quite flat; that is, they have a large radius of curvature associated with them. 
Consequently, it was difficult to determine precisely the optimum filament orientations in 
figures 10 and 11. However, it is reasonable to say from figures 10 and 11 that for large 
aspect ratios a/b > 2, 0 » ±60° to ±62° is a practical filament orientation. 

Figures 12 and 13 present interaction curves for the buckling of simply supported 
and clamped panels in combined axial compression and shear for various filament orien- 
tations and aspect ratios. The optimum filament orientations (those that correspond to 
the highest values of the buckling parameters) change according to aspect ratio a/b and 
the ratio of Nxy/Nx. For simply supported panels (fig. 12), when a/b = 1, the optimum 
orientation for all combinations of Nx and Nxy is 9 - ±45°. When a/b = 2 or 5, 
the optimum filament orientation for predominantly shear loading is near ±60° and for 
predominantly compressive loading is near ±450. For clamped panels (fig. 13) when 
a/b = 1 the optimum orientation changes from 9 = ±45° for shear loading to 9 = Qo 
for compression. When a/b = 2 or 5, the optimum orientation changes from 9 = ±60o 
for pure shear to 9 = ±45° for pure compression. This behavior was the same as that 
exhibited by simply supported panels. 
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A summary of the data from figures 12 and 13 is shown in figure 14, which indi- 
cates the banded region in which all the results lie. For orthotropic panels it was found 
that the band is bounded from below by the following simple relationship given in refer- 
ence 15 for isotropic panels: 

Rjj + Rxy (1®) 

where 


(19) 


In equations (19), Nx^ and Nxy^ are the buckling loads for pure longitudinal compres- 
sion and pure shear, respectively. Consequently, for the orthotropic cases considered, 
equation (18) is a reasonable conservative approximation for combined longitudinal com- 
pression and shear buckling of composite panels. 

Figures 15 and 16 contain, respectively, compression and shear-buckling results for 
graphite -epoxy sandwich panels with nondeflecting edge supports and r otational edge 
springs for various filament orientations and aspect ratios. The associated boundary 
conditions are given by equations (A20) to (A22), and the rotational springs were assumed 
to be uniformly distributed about the panel edges. When the spring stiffness is zero, all 
four edges are simply supported and, when infinite, all four edges are clamped. 

In general, the figures indicate that the buckling load increases sharply as the 
spring stiffness parameter bkji/Ej^t^ increases from zero to one, the buckling loads 
obtaining at least 80 percent of their clamped value when the spring stiffness parameter 
is one. With further increase in the spring stiffness the buckling loads slowly approach 
the clamped value, increasing to within at least 10 percent of the clamped value when the 
spring stiffness parameter is three. Furthermore, the curves for the ±45° filament 
orientation generally approached the clamped values most rapidly. 

CONCLUDING REMARKS. 

A computerized analysis has been developed for the combined compression and 
shear buckling of stiffened orthotropic composite panels on discrete springs. Boundary 
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conditions are general and include elastic boundary restraints. Buckling solutions are 
obtained by using a newly developed trigonometric finite -difference procedure which 
increases the solution convergence rate over conventional finite -difference methods, thus 
allowing problems to be solved with the same accuracy as with conventional differences 
but with fewer degrees of freedom. The trigonometric finite -difference procedure intro- 
duces two new parameters into the solution. These parameters can be internally selected 
by the program during problem execution or can be selected by the user. The validity of 
the program has been substantiated by comparisons with many existing known solutions. 

A program listing, input description, and sample problem are provided. 

Using the program, the classical general shear-buckling results (in terms of uni- 
versal orthotropic parameters), which are available only for simply supported panels 
over a limited rai^e of orthotropic properties, have been extended to the complete range 
of these properties for simply supported panels and clamped panels. Results for the 
shear buckling of isotropic panels with a central stiffener have also been obtained. 

The program has been applied to parametric studies which examine the effect of 
filament orientation upon the buckling of graphite -epoxy sandwich panels. From these 
studies optimum filament orientations (those which yield maximum buckling loads) were 
determined within a class of graphite -epoxy sandwich panels for all aspect ratios. In 
particular, it was found that for shear buckling of high-aspect-ratio panels (greater than 
two) reasonable filament orientations are between ±60^ and ±62° while, for axial compres- 
sion of panels with aspect ratio greater than one, a reasonable filament orientation is 
±45°. In addition, interaction curves were determined for the combined axial compres- 
sion and shear buckling of panels with varying filament orientations. A parabolic inter- 
action relationship previously developed for isotropic infinite strips in combined axial 
compression and shear provided a reasonably accurate and conservative estimate for the 
buckling loads of the orthotropic panels considered herein. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
August 1, 1975 
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DEVELOPMENT OF GOVERNING EQUATIONS 
For completeness, equations (1) to (3) of the main text are repeated here: 




(Al) 


Mx = Diiw^xx + Dl2W,yy'" 

My = Digw^xx + I^22W,yyl 

Mxy = 2D66W^xy 


(A2) 


pa pb 

= Jq Jq (NxW,x5w^x + NyW^y6w^y + NxyW^y6w^x + NxyW^x^w^y) dy dx 


(A3) 


Then, replacing the derivatives in equations (A2) by trigonometric central differences 
yields 


(w,xx)ij = AKi,) - 2wij + "1-1, j) 


^y 




y Ax Ay 


(A4) 


J 


where Ax and Ay are the trigonometric difference coefficients defined by equa- 
tion (10). The terms (w,xx)y (^,yy)y are defined at the full stations denoted by 

the circles in figure 2(b), while (w^xy) .. defined at the half stations denoted by the 
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squares in figure 2 (b). Consequently, the indices (i,j) attached to a variable may refer to 
the variable being evaluated at either full or half stations, depending on the variable. 

Introducing equations (A 2 ) and (A 4 ) into equation (Al) and replacing the double inte- 
gral by a double sum yields 

^ - 26 wij + Swi.ij) + ^ My («Wl,j+l 

Ax Ay J 


- 26 wy + 


where N and M are the total number of finite -difference stations in the x- and 
y-direetions, respectively, and Vxv ^Vi have the following definitions: 

1 I -/j 


■i) + K+i j+i - 


(A 5 ) 


6 U - Ax Ay ^ ^ 

3=1 i=l 


fo 


= \ 1/2 


vl 


(i < Ii or i > I3)' 


(i = Il or i = l3)\ 


dl < i < I3) 




(A6) 


fo 


^v. = { 1/2 


(j<J4 or i>J2)" 
(j = J4 or j = J2) y 

(J4 < 5 < J2).. 


(A7) 


fo 






(i < II or i = I3) 

(II g i £ I 3 ) 




(A 8 ) 
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(j < J4 or j £ J2) 
(J4^j<J2) 


(A9) 


In equations (A6) to (A9), Ij^ and I3 are the row designations of boundaries (D 
and (D, respectively, and 3 2 and J4 are the column designations of boundaries (D 
and 0, respectively. (See fig. 2(a).) 

Replacing the derivatives in equation (A3) by central trigonometric differences and 
the double integral by a double sum yields 


6 Vn = - 


r V r N 

= -Ax Ay 2 Z - ®"ij) * «xi’7yj 

i=l i=l ' Ax Ay 


- '>'lj)(5Wi,J + l - «W,j) + I7x,>)y, 


N 


'xy n 


3 4Ax A 


(Wi-Hlj - Wij - Wij+l)(6Wij+l 


- 6wy + 6wi+i,j+i - + (wij^.1 - Wy -t- - Wi+l,j) (6wi+i,3 - 


+ 5 wi+ij+i - 6wi,j+ij 


(AlO) 


In deriving equation (AlO), the first and second terms in the integrand of equation (A3) 
have been replaced by trigonometric differences evaluated at stations indicated by "x” 
and "y>” respectively, in figure 2(b), while the third and fourth terms have been evaluated 
at half stations, indicated by squares in figure 2(b), by averaging the derivatives. 

The external forces and moments on the panel are those comli^ from discrete 
lateral deflection and rotational springs. The virtual work of these forces and moments 
may be expressed as 

M N M N 

^ ^ kj^ijWijSwij + ^ - wij)(6Wi+ij - 6wy) 

i=lj=l i=lj=lAx 


i=l 3 = 1 Ay 


- Wij)(5wi,3+1 - Swij) 


(All) 
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where kj{i is the spring stiffness associated with a lateral deflection spring and kx 
and ky are stiffnesses associated with rotational springs which resist moments acting 
about the Y- and X-axes, respectivety. The k^ type springs act at full stations, 
indicated by circles in figure 2(b), while the kx and ky type springs act at positions 
indicated by 'Wand "y," respectively, in figure 2(b). 

Substituting equations (A5), (AlO), and (All) into the statement of the principle of 
virtual work, equation (4) yields 


M N 


Z E + %’^ij + 2NxyXjj)6wij - 0 

i=l j=l 


where 


Kii = i^x- 1 . - 2^y.Mx.. + U. ,Mx. , Mv. . , - 2^„ Mv.. 

^yj ^ 2\=xi+i sxi Xy ^ 2("’yj+i yi,j+i ^yj yy 




-%i»yj.lMxy,j.j + -7xiXyjMxy,j 


ij = ^[^yj^i^i+l,) - *ij) - «yj’’xi.j('>'li - 


^ 2 ^xi’'yj(*l,j+l - ®ij) - 5x,’)yj - '"l,j -l) 


’‘‘j " ■ "‘j) Vyj ■ ("l+l.j-1 - ®ij)’'xi’>yj.i 


- (Wij - + (wij - 
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Sij = 


k/|..W54 + 
Ax Ay ^ 


Ax Ay Ax 


‘'J'l-i.iC'ij - ’"1-1, j) - S(*«i,i ■ "ij) 


Ax Ay Ay 


“lyy-iN -SjKi+1 '"ij) 


(A17) 


From equations (A2) and (A4), the moments are related to the displacements as 
follows: 




(Mx)ij = (Dn)y(Wi^.ij - 2wij + Wl-l,i) A -1 (Dl2)ij(’»U+l - 2*ij + w,,) _i)^ 


(My) = (D22),j(wi,j+1 - 2w,j +w,_j.i)^+ (Di 2) (WH.1J - 2wy + 

A.. Ax 


/ (A18) 


(Mxy)y = 2(D66)ij(wi+i,j+i - Wij+i - Wi+i,j + wy)^ 


Ax Ay 


where (Mx).. and (Myj., act at the full stations, indicated by circles in figure 2(b), 

ij 

and (Mxy).. acts at the half stations, indicated by squares in figure 2(b). 


Boundary Conditions 


All four boundaries free or spring -supported. - If on the plate boundaries no con- 
straints exist on w or its derivatives normal to the boundary, equation (A12) must be 
valid for all virtual displacements ^wy, thus yielding equation (5) which is repeated 
here: 


^ij + ®ij + Nx«ij + ^y’^ij ■*" ^NxyXy - 0 



(A19) 


Equation (A 19) represents equilibrium at each finite-difference station with each equi- 
libriiim equation containing an array of 13 values of w as depicted in figure 2(b). In 
solving these equations by the procedure discussed in appendix C, the terms wy repre- 
sent the unknowns and equations (A18) are used to determine the moments appearing in the 
relationship for Ky, equation (A 13). 
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When a difference station lies on the boundary of the plate (that is, I = Ij or 
i = I 3 or 3 = J 4 or 3 = J 2 ), the corresponding equilibrium equation reduces to the nat- 
ural boundary condition on the Kirchhoff shear, reference 6 . Also, when a difference 
station lies one finite difference interval off the plate (that is, i = Ij - 1 or I = I 3 + 1 
or 3 = J 4 - 1 or 3 = J 2 + 1 ), the corresponding equilibrium equation reduces to the nat- 
ural boundary condition on the bending moment. Furthermore, when a difference station 
lies two or more finite-difference intervals off the plate (that is, i < I 4 - 1 or i > I 3 + 1 
or 3 < J 4 - 1 or 3' > J 2 + 1 ), the corresponding equilibrium equations reduce to the triv- 
ial equation 0 = 0. Consequently, no equilibrium equations exist for these stations. 

Edges with nondeflecting lateral supports and rotational springs . - Equation (A19) 
may be used in approximating the solution of problems with nondeflecting edges; for 
example, if w = 0 on an edge, equation (A19) may be used in conjunction with extremely 
stiff lateral springs placed along the edge. Alternatively, an edge which is restrained 
from lateral motion may be handled as a special case, and in so doing the number of com- 
putations required for the problem solution is reduced. 

The boundary condition for a nondeflecting edge is 

w = 0 (on the edge) (A20) 

If, in addition, uniformly distributed rotational springs act along boundaries (I) and (D 
(see fig. 2 (a)), 

Mx = kRW^x (on the edge) (A21) 

or, if uniformly distributed rotational springs act along boundary or 

My = kRW^y (on the edge) (A22) 


As a result of the foregoing, equation (A20) replaces the boundary condition on the 
Kirchhoff shear, while the difference form of equation (A 21 ) or (A22) replaces the bound- 
ary condition on the edge moment. Furthermore, as an example, equation (A21) on 
boundary (T) becomes 





(A23) 
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where 


A * TTAy 

2\xSin^-^ 


(A24) 


Substituting for Mx from equations (A18) and employing equation (A23) yields 


Then 


V i K.R 


(A25) 




Al“ a,2 


w 


Il+l,j 


kt. i 

*^R 

>1 * o 

Ax Ax 


(A26) 


Substituting into the first of equation (A25) yields 








1 + 



(Dll)l, 

% 

■ l^j 

Ax 

1 2 

Ax 




A,2 _ 


w- 


il+1,3 


(A27) 


It is evident from an examination of the first of equation (A25) that equation (A23) is 
satisfied by setting = 0 and (Dn)j^ . = . where 




1 + 




A^ 

Ax^ 

^ + 


A* 

Ax^ 


(Dll) 


Il>3 


(A28) 


Similar relationships maybe developed for boundaries (|), d), and 0, 


23 



APPENDIX A 


In summary, for a nondeflecting boundary with uniformly distributed rotational 
springs, equilibrium on the boundary and one station off the boundary are not used. 
Instead, in the remaining equilibrium equations, w on the boundary and one station off 
the boundary are set equal to zero and on the boundary is set equal to D|j^ if the 

boundary is number (J) or (D, and D 22 on the boundary is set equal to D |2 if the 
boundary is number d) or 0 , 

The limiting cases of simply supported or clamped boundaries are readily provided 
by letting kj^ approach zero or infinity, respectively. Hence, for a simply supported 
boundary 

= 0 if the boundary is (J) or (|) 

D |2 = 0 if the boundary is or @ 
and for a clamped boundary 


= 2Djj if the boundary is (|) or (D 
D ^2 = 2 D 22 if the boundary is (D or 0 


Flexural Stiffeners 

The effects of flexural stiffeners are accounted for in a manner similar to that used 
for nondeflecting supports. At each finite -difference station along the stiffener, (Dll)y 
is replaced by (Du).j if the stiffener is parallel to the X-axis and (D22)ij replaced 
by (D22)j^j if the stiffener is parallel to the Y-axis, where 


> 

(D22)y = (D22),j * S 


(A29) 


and El is the lateral bending stiffness of the stiffener about the neutral plane of the 
panel. 
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Summary of Finite -Difference Stations at Which Equilibrium Is Enforced 

As a result of the foregoing discussions on free or spring -supported edges and non- 
deflecting e<%es, the rows i and columns j at which equilibrixxm is enforced are, 
respectively, 


(A30) 


Ix + 3 - Twice the number of nondeflecting edges parallel to the Y-axis' 


Ng = J 2 - J 4 + 3 - Twice the number of nondeflecting edges parallel to the X-axis 
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TRIGONOMETRIC FINITE DIFFERENCES 

Trigonometric finite differences introduce the trigonometric parameters 
\y which are not present in conventional finite differences. Consequently, the first pur- 
pose of this appendix is to present and demonstrate some effective procedures for 
selecting values of A-x \y which results in an improved convergence rate over 

conventional differences. The second purpose is to point out some of the limitations of 
trigonometric finite differences. 


Selection of and \y 

Selection of values of A-x a^nd Ay which improve the convergence rate of trigo- 
nometric finite -difference solutions over those of conventional finite -difference solutions 
is predominantly based on engineering considerations and experience. Experience has 
shown that it is often advantageous to select trigonometric parameters whose ratio is 
determined on the basis of the infinitely long panel solution as is done in equations (11) 
and (12), that is. 



(Bl) 


where /3 is the wavelength parameter of an infinitely long panel, defined as the ratio 
of the panel width to the buckle lei^th. Imposing equation (Bl) on the parameter selec- 
tion should be reasonable for panels which buckle with more than two half waves along 
their length. 

The value of /3 may be determined to any degree of accuracy by extending the 
isotropic results of reference 16. For a panel with its long dimension parallel to the 
X-axis, first approximations of the buckling eigenvalue Poo and wavelength parameter 
j3 satisfy the following two simultaneous equations for panels whose long sides are sim- 
ply supported: 


'\ 


(txy + Poo^^xy)^ - I M 1 M 2 = 0 
^(MlMg) = 0 




(B2) 
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and, for panels whose long sides are clamped, p^ and 13 satisfy the two simultaneous 
equations 


(txy + P„rxy)2 - ||(2Mo + M 2 )(Mi + Mj) = 0 
i(Mo + M2 )(Mi + Mj) = 0 


(B3) 


where 




^ ^ + Poo^x) - n2(ty + p^ry) 


Dll 


(n = 0, 1, 2, 3) 


Convergence Behavior 

Figures 17(a) to 17(f) illustrate the convergence of trigonometric finite -difference 
solutions when is fixed on the basis of equation (Bl). Results for both simply 

supported and clamped isotropic panels under either axial compression or shear are 
shown in these figures. In each case the panel was modeled using an equal number of 
finite -difference stations in the x- and y-directions. Exact and approximate values for 
these cases are given in references 1, 6, 16, and 17. 

The dashed curve in each of figures 17(a) to 17(f) illustrates the convergence of the 
conventional difference solution - that is. Ax and Ay infinite - while the solid and 
dash -dot curves illustrate the convergence achieved with some finite values of ^x. Com- 
parison of the curves indicates that some values of Ax increase the convergence rate 
over the conventional rate while other values decrease it. (In those special cases where 
the buckle shape is exactly a double sine wave, the trigonometric difference solution is 
exact when Ax and Ay are equal to the buckle half wavelength.) Consider though the 
dash-dot curve of each figure. These curves show the convergence when Ay is simply 
taken equal to the panel width and Ax is taken equal to the buckle length of the infinitely 
long panel; that is, equations (11) and (12) are applied. Comparison of the dash-dot 
curves and the dashed curves indicates that equations (11) and (12) provide reasonable 
values of Ax and Ay which improve the solution convergence. As figures 17(a) to 17(f) 
indicate, however, other values of Ax /a could be selected which further improve the 
convergence rate. Such values may be found by making a condensed cross plot of each 
figure; for example, consider the case of the compression of a square isotropic clamped 
panel as shown in figure 17(c). For this case, equations (B3) predict 0= 1.5. Then, 
using program BOP with Ayy/Ax = 1.5, Ax /a is varied from 0.25 to 1 for mesh sizes of 
a/Ax = b/Ay = 5 and a/Ax = b/Ay = 6; these curves are shown in figure 18. As the 
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mesh spaeir^ is decreased, the curves will approach the exact solution at all values of 
A-x/a. However, the two curves cross at Ax/a = 0.35 and Nx = 9.75, which implies 
that convergence is most rapid at this value of A-x/a since increasing the mesh size 
did not change the buckling stress coefficient. It is evident from figure 17(c) that, if 
such a choice of A.x were used, convergence would be improved beyond that achieved 
by selecting A.^ from equation (11), 

As further examples, consider the results in table 4 for the shear buckling of the 
orthotropic panels described in table 3. The values of A.x and A.y were determined 
by making the required cross plots. It is evident by comparing the conventional and 
trigonometric solutions given in the table that the selected values of A,x and A.y pro- 
vided excellent results. 

The additional effort involved in finding better values of Xx be justified in 

problems where convergence would otherwise be extremely slow. It may also be justi- 
fied in the performance of parameter studies. In such studies some t 3 rpical problems 
within the problem class to be studied are chosen; for these, improved values of Xx 
found and then interpolated to yield Xx for other problems within the study class. 


Correction Factors for Equations (11) and (12) 

Equations (B2) and (B3) which provide jS for equations (11) and (12) do not cover 
every case; the boundary conditions may not be simply supported or clamped, or it may 
be inappropriate to use |8 based on an infinitely long panel. Consequently, equa- 
tions (11) and (12) must be used with engineering judgment. Some allowance is provided 
by introducing correction factors Cx and Cyx into equations (11) and (12), that is. 



(B4) 


^x _ b ^x 
a a /3 


(B5) 


A numerical routine which calculates ^ from equations (B2) or (B3), and then ^x 
and Xy from equations (11) and (12), is used in program BOP. This program is briefly 
discussed in the main text and is documented in appendix D. 


Limitations of Trigonometric Finite Differences 

In figure 19 a sketch of the variation with Xx of the coefficient i/Ax as defined 
by eqmtion (10) is presented. The reader's attention is called to the singularities of 


28 



APPENDIX B 


1/Ax at \x = 6tc. In order to avoid these singularities and the rapidly 

Z 4 o 

varying behavior of i/Ax between them, \x 3^nd similarly Xy must be chosen such 
that 




Xx > 




(B6) 


Moreover, if uniformly distributed rotational springs are prescribed on the boundaries in 
the manner presented in equations (A20) to (A24), then to avoid singularities in A^ and 
A^ choose 


^ ^x 


Xy > Ay 


(B7) 
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STABILITY DETERMINANT EVALUATION 

Since the total number of rows and columns at which equilibrium is enforced is 
Me and Ne, respectively, a stability determinant of order MeNe x MeNe would result. 
To produce a stability determinant of smaller size, a marching procedure is employed. 
This procedure, which is described herein, operates on the equilibrium equations to pro- 
duce, by a process of successive elimination, a determinant of size 2Me x 2Mg. 

The marching procedure takes advantage of the fact that each of the difference 
equations of equilibrium, equations (5), is linear arid homogeneous, with each one con- 
taining no more than 13 unknown deflections. For a station (i,j) away from the plate 
edges 


If + l<i<l£-l 
+ -1 

where If and Ij^ are the first and last rows of finite-difference stations at which 
equilibrium is prescribed, and Jf and Jj^ are the first and last columns of finite - 
difference stations at which equilibrium is prescribed, the 13 unknown deflections form 
the geometric pattern shown in figure 2(b). It is evident from this pattern that the deflec- 
tions at stations in column j + 2 can be determined by using equilibrium at stations in 
column j if the deflections in columns j - 2, j - 1, j, and j + 1 are known or pre- 
scribed. For equilibrium at stations lying near the edges, however, the geometric pat- 
tern of figure 2(b) is reduced. Consequently, equilibrium at stations in the first column 
Jf may be used to determine the deflections at stations in column Jf + 2 if the deflec- 
tions only in columns Jf and Jf + 1 are prescribed, since deflections in columns 
Jf - 1 and Jf - 2 do not appear in these equilibrium equations. 

Having found the deflections in column Jf + 2 from prescribed values in column 
Jf and Jf + 1, equilibrium at stations in column Jf + 1 can be used to obtain the deflec- 
tions in column Jf + 3; likewise, equilibrium at stations in column Jf + 2 can provide 
deflections in coliunn Jf + 4, etc. Thus, a marching routine is developed from column 
to column which determines the deflections throughout the panel from prescribed values 
in the first two columns. It should be noted that equilibrium at stations in the last two 
columns, J^ - 1 and Jj^, is not used at this stage of the marching procedure. 

The evaluation of the stability determinant can now be performed numerically for a 
given value of the eigenvalue by choosing 2Me linearly independent sets of assumed 
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deflections for the first two columns. 


These assumed sets are taken as 



where each column contains 2M© values. By marching across the plate with the rth set 
of these assumed values, deflections throughout the plate wj?^ are determined. How- 
ever, the equilibrium equation at stations in the last two columns will not, in general, be 
satisfied by any of these assumed sets. Therefore, consider the column matrix 


e 


(r) 

If,j£-1 



(C2) 


h’h 


where each element of the matrix represents the value, of the left-hand side of an equilib- 
rium equation at a station in columns - 1 or for the rth assumed set and would 
be identically zero if the assumed deflections were exact. The total solution is a linear 
superposition of all the assumed sets, that is. 



(C3) 


Jf ^ j g J£ 
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Correspondingly, the total contribution to equilibrium at columns - 1 and Jjg for 
all assumed sets of deflections is 

2Me 

[e] = ^ (C4) 

The coefficients A(r) 

are determined by enforcing equilibrium at stations in the last 
two columns which leads to 

[e] = 0 (C5) 



For a nontrivial solution of equation (C6) the determinant of the coefficients must vanish, 
resulting in 

|e| ^0 (C7) 

and it is clear from equation (C6) that |e| is of order 2Me x 2Me. 
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COMPUTER PROGRAM 

The computer program BOP (Buckling of Orthotropic Panels) was written in 
FORTRAN IV on a SCOPE 3.1 system modified for Langley Research Center and executes 
and loads with a field length of 60000 octal locations. The program is applicable to the 
combined compression and shear of stiffened, variable -thickness, flat rectangular ortho- 
tropic panels on discrete springs; boundary conditions are general and include elastic 
boundary restraints. A description of the input, an example problem showing input and 
output, and a program listing are provided. 

Input Description 

For each case the input consists of a single identification card and a Namelist 
BUCKLE as follows: 

ISTIFF,ISTEP,IX,JX,MSHAPE,MA,NOMAT,TH,AT,MATYPE,El,E2,Ul,G12,IBC,AKR,Dl, 
D2,D12,D66,DS1,XA,XB,AKL,AKX,AKY,NUPRIT,EI,I0RIENT,LCX:,TX,TY,TXY,RX,RY, 
RXY, PI, DELP,PFIN,TEST,MR,NC,X, Y, DS2, DS 12, DS66 

Many of the input variables have associated default values as will be indicated in the fol- 
lowing descriptions: 

Control parameters 

ISTIFF =1 no preprocessing of laminate properties - execute for buckling (only) 

= 2 preprocess and execute for buckling 
= 3 preprocess only — do not execute for buckling 
DEFAULT :ISTIFF = 2 

ISTEP =1 prc^ram automatically varies the input step size, DELP 
= 2 step size fixed and equal to DELP 
DEFAULT:ISTEP = 1 

tx = 1 output of intermediate results 

= 2 oulput of intermediate results suppressed 

JX = 1 output of flexural stiffnesses at each finite -difference station 

= 2 output of flexural stiffnesses suppressed 

DEFAULT:IX = JX = 2 
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MSHAPE = 1 compute mode shape 

= 2 do not compute mode shape 

DEFAULT:MSHAPE = 2 

Laminate and lamina properties (Required if ISTIFF = 2 or 3) 

MA number of laminas in the laminate 

NOMAT number of different materials comprising the laminate 

TH a one -dimensional array in which the ith element of the array corre- 

sponds to the filament orientation (as measured from the X-axis in 
degrees) in the ith lamina 

AT a one -dimensional array in which the ith element of the array corre- 

sponds to the thickness of the ith lamina 

MATYPE a one -dimensional array in which the ith element is the number desig- 
nation of the material in the ith lamina 

El a one -dimensional array in which the jth element of the array corre- 

sponds to the Young's modules parallel to the fibers in the jth 
material 

E2 a one -dimensional array specifying the Young's modulus transverse to 

the fibers 

U1 a one -dimensional array specifying Poisson's ratio in each 

lamina 

G12 a one -dimensional array specifying the shear modulus in each material 

Boundary conditions 

IBC a one -dimensional array of four elements in which the ith element 

refers to the ith boundary (see fig. 2(a)); four options are available 
at each boundary 

IBC(I) = 1 nondeflecting lateral support with uniform rotational springs on edge I 
= 2 simple support on edge I 
= 3 clamped on edge I 
= 4 free on edge I 

= 5 other boundary conditions - set by user through appropriate input of 
Dl, D2, D12, and D66 
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AKR a one -dimensional array in which the ith element of the array corre - 

spends to the uniformly distributed rotational spring stiffness per 
unit length of boundary on the ith boundary; required if any boundary 
has IBC = 1 

Laminate flexural stiffnesses (Required if ISTIFF = 1) 

D1 a two-dimensional array in which the (i,j)th element of the array cor- 
responds to the value of y 

D2 similar to Dl, but specifying (D22)ij 

D12 similar to Dl, but specifying (Di2)y 

D66 similar to Dl, but specifying (D66)ij 

DSl reference value of Dj^j^ 

Plate geometry 

XA = a dimension parallel to X-axis (fig. 2(a)) 

XB = b dimension parallel to Y-axis (fig. 2(a)) 

Discrete springs 

AKL a two-dimensional array in which the (i,j)th element corresponds to 

AKX similar to AKL but referring to (hx)^j 

AKY similar to AKL but referring to _ 

Discrete flexur al stiffeners 
NUPRIT 

El 

lORIENT 

= 1 
= 2 

LOG 

DEFAULT :NUPRIT = 0; El, lORIENT and L(X! need not be input 


number of stiffeners 

a one -dimensional array whose ith element specifies the flexural stiff- 
ness of the ith stiffener about the neutral plane of the panel 

a one -dimensional array whose ith element specifies whether the stiff - 
ener is parallel to X- or Y-axis 

stiffener parallel to X-axis 

stiffener parallel to Y-axis 

a one -dimensional array whose ith element gives the row or column 
location of the ith stiffener 
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Applied in -plane loads 

In -plane loads are assumed to be uniform over the boundary to which they are 
applied and are increased to buckling according to the relationships prescribed by equa- 
tions (13); therefore, the user inputs 

TX = Fx 

TY=Fy 

TXY=Fxy 

RX=rx 

RY=ry 

RXY = rxy 

Eigenvalue search parameters 

PI starting value of p. If PI < 0., the program will calculate PI from 

equation (B2) or (B3) according to the relation, 

PI = ABS(P 1) *PBAR (Dl) 

where PBAR is p^ from equation (B2) or (B3) , 

DEFAULTrPl = 0.9*PBAR 

DELP increment of (p); if PI < 0., DELP = 0.1+PBAR; if ISTEP = 1, 

DELP is automatically varied during the eigenvalue search 

PFIN maximum value of p during the eigenvalue search 

TEST eigenvalue accuracy 

DEFAULT:!, x 10 -3 

Trigonometric f inite -difference data 

MR number of rows of finite -difference stations interior to the plate ^ not 

including boundaries 

NC number of columns of finite -difference stations interior to the plate — 

not including boundaries 

Note: The marching procedure requires NC & 4 
X = ^-x/a 

Y = \y/b 
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Note: If the user inputs X g 0, the program automatically calculates a new value of X 
and Y according to the relationship expressed by equations (B4) and (B5); that is, 

X = ABS(X) *XB/BETA/XA (D2) 

Y = ABS(Y) (D3) 

where the input magnitudes of X and Y (that is, ABS(X) and ABS(Y)) replace 
and Cyx in equations (B4) and (B5). Also, in equation (Dl), BETA = /3, and /3 is cal- 
culated from equation (B2) or (B3). 

When ISTIFF = 1 and the evaluation of X and Y is chosen, the user must also 

input 

DS2 average or typical value of D 22 

DS12 average or typical value of X>i2 

DS 66 average or typical value of Dg 0 

DEFAULT: Calculation of X and Y using equations (D2) and (D3) where ABS(X) 
and ABS(Y) are set equal to unity. 

Example Problem 

Consider the shear buckling of a 12 -inch by 3 -inch clamped sandwich panel which 
has as its lay-up, 45/-45/45/-45/CORE/-45/45/-45/45. The core thickness is 
0.0605 inch and each lamina of the skins is graphite -epoxy with a thickness of 
0.0055 inch. 

Sample Input 

THfb IS A i-'RcE FIELD I OEnT 1 F I QATl ON CARO 
? IS U C K L £ T X = . 0 , T Y = . 0 , T X Y = . 0 , R X - . 0 , R Y = . 0 , R X Y = 1 . , 

XA=1 2 , Xii = i . /MR = 1 2 , NC=6 , I 6C=R*5 , NuMA T = 2 , E 1 = 2 , 1 0 E 7 , 1 . , £ 2 = 2 . 59 £6 /I . , 

Ui = , 31, .2,G12=6.5t5, 1 . ,MA=i),MATyPE = 4*l , 2 , 4*1, 

AT = 4* .U035, ,060 '3,4* . 0035, Tti»=4:3 . ,-45. , 43, , -4 5 . , . 0, -4 5 . , 45 . , -4 5 , , 45 , 
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Sample Output 


INPUT FOR CASE 


THIS IS A FREE FIELD INDENT I FICAT ION CARD 


4 2 

4 2 

4 2 

4 2 

4 2 

4 2 

333333333333 

A 

★ 

* 

* 


X 


■ ■ISTIFF = 2 TSTEP = 1 ’ IX = 2 JX=2 TEST= 

X^A=--‘T 3.00000000F>00 ASPECT 

MR=12 NC= 6 

TX= 0. TY= 0. TXY= 

Pl= -9. OOOOOOOOE-01 DELP=" I .OOOOOOOOE-01 PFIN 

x="’^iraaooooooEFoo y= i.ooooooooE^'Oa 

RX^ 0. RY= 0. RXY= 


l.OQOOOOODE-03 

RAT 10= 4400000000E4'00 

0 . 

^ leOOOOOOOOE+02 
l.OOOOOOOOE+00 
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LAMINATED PLATE PROPERTIES 


MATERIAL KIND El E2 U1 GXY 

1 2. lOOOOOOOE <-07 2. 39000000E+ 06 3 . lOOOOOOOE-01 6.50000000E4-05 

■ - - 2 l.OOOOOOOOE+DO 1. OOOOOOOOE<-O0 2 .OOOOOOOOE-01 I . OOOOOOOOE+00 

“TAV'EfT'TjO MAT. KIND THICK THETA ' 

1 1 5 .50000000E-03 4.50000000E4-01 

2 1 5.5000000DE-03 -4. 50030000E+01 

3 I ' 5.50D00D00F-03 4.50000000E+01 

4 1 5.50000000E-03 -4. 50000000E+01 

5 2 6.05 000000E'-02 0. 

6 “ - • - -X ■5.50000000E'-03 -4. 30000000E +01 

7 1 5.5OOOOOO0E-O3 4. 50000000E +01 

8 1 5.50000000E-03 -4. 50000000E +01 

“g -i' - 5.50DDO000F-O3 4.500‘OOOOOF+OI 


^ 1 **+#***#+***+ ***********+*****#**#* **♦****#♦*♦*#**#****#***♦****♦**** ****** 


CAUTION — COUPLING BETWEEN EXTENSION AND BENDING MAY BE SIGNIFICANT 
IF THE FOLLOWING FOUR VALUES ARE NOT ALL EQUAL. 

XF THIS IS' THE C AS E V THE" RESULTS SHOULD' BE USED WITH DESCRETION ' 

4.76775528E-17 7. 15163293E-17 8. 80l022g6E-lT 5.97234077E-17 

35c*jjc* ****** :5c )Sc;}!:5c79c7St 


A MATRIX 


3,052152295 + 05 
-:2 . W0r51T9E+XJ5 

0 , 


2. 480151 79E +05 
3.052152 295+1)5 

0 . 


0 , 

0 . 

2.436551 39E +05 


B MATRIX 


1.45519152E-11 
2'.T82T8T25F~11 
7. 2 759576 IE- 11 


2. 182787285-11 
2.18278T2S5--1I' 
7. 275957615-11 


8.003553385-11 

T.2T5957615-11 

1.45519152E-11 
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0 MATRIX 


5.31653046E + 02 4. 320 16589E 4-02 4.69571561E+01 

4.3201 658 9E 4-02 ’ 5. 316 530 46E4-02 ' 4.'6957 1 56 lE+0 1 

4. 69571561E4-01 4. 69571561E4-01 4.24421785E4-02 


LAMI];^AJE properties 

EX= 9.92156003E4-05 EY= 9.92156003E4-05 GXY~ 2.33162813E4-06 
NUXY= 8.12591100E-01 NUYX= 8. 12 5 911 OOE -0l 


P= 1.30776204E4-01 B= 9.57954000E-01 F- -1. 83044904E-08G= 3 . 83 577069E-03 


PROGRAM HAS COMPUTED AND USED X= 2.60972865E-01 

AND' Ys" 'X;'0D'D0'00'00E+DD ■8ASE0 THE INFINITE PLATE WAVE LENGTH 


■ BASE'D W'AN 'IWINr TEL Y LONG plate THE euCKLlNG STRESS COEFFICIENTS 'ARE 
NX6AR= 0. 

NY8AR= 0. 

NX7ffA'R= “IVIOTT^^ 

AND THF STRAINS ARE 
STRNX= 0. 

STRNY= 0." "" 

STRNXY= 3. 12923370E-02 

PROG'RXM“ WTLr'N0'W"'CONTTN FINITE ASPECT RATTO SOLUTION ' ' ' 


♦ :Jc:<c3{c **5!c * jjt 5}e** 


BOUNDARY CONDITIONS 
BOUNDARY NO. 1 IS CLAMPED 


BOUNDARY NO. 2 IS CLAMPED 


BOUNDAPY NO. 3 IS CLAMPED 
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BOUNDARY NO. 4 


IS CLAMPED 


T^XB AR ■ 


NY BAR 


NXYBAR 


DETERMINANT 


0 . 

■D." 

0 . 

0 . 

0 . 

0 . 


0 . 
0 . 
0 . 
0 . 
0 . 
0 . 
0 . 


1. 176 98 583 E +01 
1. 1 9 006345 E+ 01 
1.20314107E+01 
1. 20473676 E+01 
1.20687996E+01 
1. 20 761 369 E+01 
1.20793082E+01 


6.54187653E+12 
2.18161022E+rZ 
2. 37244168E+11 
1.33435 159E+11 
3.35545245E+I6 
9.37018457E+09 
0 . 


"NX= 0.~ 


BUCKLING LOADS PER UNIT OE LENGTH ALONG BOUNDARY EDGE 

NY= 3. NXY= 7.04251216E+D3 


IXB**2) *7sIX/{PI=^^*2)/{T*«3 J = 0. 

1XB**2) WY/(PI**2)/(T**3) = 0. 

'TXB^^?r'^NXT7rRT*^^2T7rT^«^ ) = ^ ‘ 5 . 6275 776 7ET06 


(XB**2 ) ^EPSILONX/ CT*=?^2) 
rXB-^^2 ) ^PSTLONT/ f ^21 
( XB**2 ) *EP SI LONXY/ (T**2 


- 0 . 

-1.19105540E+01 
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13200000 
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TABLE 1.- SHEAR-BUCKLING LOAD COEFFICIENTS FOR RECTANGULAR ORTHOTROPIC PANELS 
WITH ALL EDGES SIMPLY SUPPORTED AND THE TRIGONOMETRIC DIFFERENCE 
PARAMETERS ON WHICH THEY ARE BASED 


Stiffness parameter, Aspect -ratio parameter, 

^ \f^U^ „ b\4/^ y-directions 

a/ Ax V^y 


Wavelength ratios 
used in trigonometric 
differences 


^x/3- 

Xy/b 

.2^Di 

0.56 

0.56 

26.28 

.56 

.60 

21.43 

.56 

.80 

17.33 

.50 

.90 

15.36 

.50 

1.00 

13.77 

.35 

1.00 

11.55 



10.87 

.56 

.56 

15.78 

.56 

.60 

12.98 

.56 

.80 

10.86 

.50 

.90 

9.93 

.50 

1.00 

9.29 

.30 

1.00 

8.21 

--- 

--- 

7.72 

.56 

.56 

12.21 

.56 

.60 

10.11 

.56 

.80 

8.67 

.50 

.90 

8.09 

.50 

1.00 

7.73 

.25 

1.00 

6.71 



6.53 

.56 ' 

.56 

10.40 

.56 

.60 

8.66 

.56 

.80 

7.57 

.50 

.90 

7.10 

.-50 

1.00 

6.80 

.25 

1.00 

6.02 


— - 

5.79 

.56 

.56 

9.31 

.56 

.60 

7.68 

.56 

.80 

6.91 

.50 

1.00 

6.22 

.23 

1.00 

5.49 



5.33 


‘'For B =s 0, kg was calculated by using equations (B2). 




TABLE 1. - SHEAR-BUCKLING LOAD COEFFICIENTS FOR RECTANGULAR ORTHOTROPIC PANELS 
WITH ALL EDGES SIMPLY SUPPORTED AND THE TRIGONOMETRIC DIFFERENCE 


PARAMETERS ON WHICH THEY ARE BASED - Concluded 



^ For B = 0, ks was calculated by using equations (B2). 





TABLE 2.- SHEAR-BUCKLING LOAD COEFFICIENTS FOR RECTANGULAR ORTHOTROPIC PANELS 
WITH ALL EDGES CLAMPED AND THE TRIGONOMETRIC DIFFERENCE 
PARAMETERS ON WHICH THEY ARE BASED 


Stiffness parameter, Aspect-ratio parameter, 
JD11D22 _ _ b\4/DiT 


No. of mesh points 
in x~ and 
y -directions 


Wavelength ratios 
used in trigonometric 
differences 


Shear -buckling load 
coefficient, 


ay ^22 


b/Ay 

Xx/a 

Xy/b 


1.0 

9 

9 

1.00 

1.0 

32.56 

.8 

9 

9 

.80 

1.0 

26.31 

.6 

9 

9 

.60 

1.0 

22.21 

.4 

11 

11 

.40 

1.0 

18.91 

.2 

17 

9 

.31 

1.0 

17.34 

.1 

25 

9 

.15 

1.0 

17.31 

0 

— - 


--- 


17.13 

1.0 

0 

9 

1.10 

1.1 

21.63 

.8 

9 

9 

.90 

1.0 

17.92 

.6 

9 

9 

.60 

1.0 

15.43 

.4 

11 

11 

.40 

1.0 

13.62 

.2 

17 

9 

.25 

1.0 

12.64 

.1 

25 

9 

.13 

1.0 

12.89 

0 



--- 

--- 

12.51 

1.0 

9 

9 

1.10 

1.1 

17.86 

.8 

9 

9 

.90 

1.0 

14.89 

.6 

'9 

9 

.60 

1.0 

13.06 

.4 1 

11 1 

11 

.40 

1.0 

11.60 

.2 

15 

8 

.22 

1.0 

10.64 

.1 ^ 

25 ! 

9 

.13 

1.0 

10.95 

0 





10.69 

1.0 

9 

9 

1.10 

1.1 

15,94 

.8 

9 

9 

.90 

1.0 

13,34 

.6 

9 

9 ■ ! 

.60 

1.0 

11,84 

.4 

11 

11 

.40 

1.0 

10.55 

.2 

17 

9 

.24 

1.0 

9.99 

.1 

25 

9 

.13 

1.0 

10.16 

0 


— - 


— - 

9.63 

1.0 

9 

9 

1.20 

1.2 

14.81 

.8 

9 

9 

1.00 

1.0 

12.44 

.6 

9 

9 

.60 

1.0 

11.08 

.4 

11 

11 

.40 

1.0 

9.89 

.2 

17 

9 

,22 

1.0 

9.27 

.1 

25 

9 

.12 

1.0 

9.11 

0 




— 


8.99 


^For B = 0, ks was calculated by using equations (B3). 




TABLE 2.- SHEAR-BUCKLING LOAD COEFFICIENTS FOR RECTANGULAR ORTHOTROPIC PANELS 
WITH ALL EDGES CLAMPED AND THE TRIGONOMETRIC DIFFERENCE 
PARAMETERS ON WHICH THEY ARE BASED - Concluded 


Stiffness parameter, Aspect -ratio parameter, Inx^and^^^^^ 

\/DiiD 22 g _ b\4 /Pll y -directions 

® ” TN u\/ Doo — — r — — 


Wavelength ratios 
used in trigonometric 
differences 


Shear -buckling load 
coefficient. 



a/ Ax 

b/Ay 

^x/9- 

Xy/b 


1.0 

9 

9 

1.20 

1.2 

13.87 

.8 

9 

9 

1.00 

1.0 

11.68 

.6 

9 

9 

.60 

1.0 

10.46 

.4 

9 

9 

.40 

1.0 

9.39 

.2 

15 

8 

.22 

1.0 

8.80 

.1 

22 

8 

.12 

1.0 

8.98 

0 




--- 

8.45 

1.0 

9 

9 

1.20 

1.2 

12.91 

.8 

9 

9 

1.00 

1.0 

10.90 

.6 

9 

9 

.60 

1.0 

9.80 

.4 

9 

9 

.40 

1.0 

8.86 

.2 

15 

8 

.22 

1.0 

8.34 

.1 

22 

8 

.12 

1.0 

8.58 

0 

— - 

--- 



7.93 

1.0 

9 

9 

1.20 

1.2 

11.93 

.8 

9 

9 

1.00 

1.0 

10.11 

.6 

9 

9 

.60 

1.0 

9.07 

.4 

9 

9 

.40 

1.0 

8.31 

.2 

15 

8 

.22 

1.0 

7.84 

.1 

25 

9 

.12 

1.0 

8.12 

0 


--- 

— - 

— - 

7.32 

1.0 

9 

9 

1.20 

1.2 

10.94 

.8 

9 

9 

1.00 

1.0 

9.31 

.6 

9 

9 

.60 

1.0 

8.33 

.4 

9 

9 

.40 

1.0 

7.74 

.2 

15 

8 

.22 

1.0 

7.33 

.1 

25 

9 

.12 

1.0 

7.66 

0 

— 


— - 

-- 

6.72 

1.0 

9 

9 

1.20 

1.2 

9.92 

.8 

9 

9 

1.00 

1.0 

8.48 

.6 

9 

9 

.60 

1.0 

7.57 

.4 

11 

11 

.40 

1.0 

6.97 

.2 

15 

8 

.22 

1.0 

6.79 

.1 

25 

9 

.12 

1.0 

7.17 

0 



--- 


6.11 


^For B = 0, kg was calculated by using equations (B3). 



TABLE 3.- MATERIAL PROPERTIES OF GRAPHITE -EPOXY SKINS 
WITH THEIR EQUIVALENT ORTHOTROPIC PARAMETERS 
AT VARIOUS FILAMENT ORIENTATIONS 

I^Ei = 145 GN/m2 (21 x 106 psi); E 2 /E 1 = 0.1138; 


Gi2/Ei = 

0,03095; 1^12 = 0-31J 


Filament orientation, 
±0 , deg 

_ \/DiiD22 ^ 
D 3 

a B = 

b ' V®22 

0 

3.50 

1.722 

30 

.511 

1.389 

45 

.415 

1,000 

60 

.511 

.720 

90 

3.50 

.581 




TABLE 4.- COMPARISON OF CONVENTIONAL AND TRIGONOMETRIC FINITE DIFFERENCES 


CQ 

m 

52 ; 

<d 

ft 

u 

s 

ft 


ft 

o 

o 
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Figure 1.- Stress resultants acting upon an 
element of the plate. 
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Evaluation of 



X j-2 j-1 j j+1 

(a) Panel geometry and boundary (b) Finite -difference station layout and 

designation. designation. 



Figure 2. - Geometrical and numerical configurations. 
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Shear-buckling load coefficient, 






Shear-buckling load coefficient 





igure 
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Filament orientation, ± 6 

Figure 6.- Variation of compressive buckling load with filament orientation for panels of various aspect ratios. 



Shear buckling parameter 


b 



Filament orientation, ± 6 

Figure 7. - Variation of shear buckling load with filament orientation for panels 

of various aspect ratios. 
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Aspect ratio, a/b 

rientation for the compressive buckling of a simply supported panel. 


— 1 
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Figure 9>- Optimum filament orientation for the compressive buckling of a clamped panel. 









, deg 



Figure 13. - Combined axial compression 






Figure 14,- Summary of combined axial compression and shear-buckling 
results for simply supported and clamped panels. 
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(b) a/b = 2. 
Figure 15 .- Continued. 





± e , deg 



Figure 16.- Variation of shear buckling parameter with rotational edge spring stiffness. 




± 6 f deg 



Figure 16. - Continued. 





Figure 16*- Concluded. 




Percent 

error 



0 20 40 60 80 100 


Degrees of freedom 


(a) Convergence for the eompressive buckling of simply supported isotropic 
square panels, ^y/^x = ^ = 1* From reference 6, Nx = 4.0. 

Figure 17.- Convergence characteristics of trigonometric finite differences. 
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Degrees of freedom 


(b) Convergence for the compressive buckling of 5 X 1 simply supported isotropic panels. 

^y/^x = ^ = 5. From reference 6, Nx = 4.0. 

Figure 17. - Continued. 
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Percent 

error 


0 


.33 



0 20 40 60 80 100 
Degrees of freedom 

(c) Convergence for the compressive buckling of clamped square isotropic panels 
\y^\x = = 1.5. From reference 17, Nx ~ 10.074. 

Figure 17.- Continued. 





(d) Convergence for the compressive buckling of isotropic clamped 5x1 panels, 
^y/^x = 0 From reference 16, Nx * 7.0. 

Figure 17.- Continued. 




Conventional 



0 40 80 120 160 200 
Degrees of freedom 


(e) Convergence for the shear buckling of isotropic singly supported 5X1 panels 
Xy/Xx = ^ = 0.8. From reference 1, Nxy “ 5.55. 

Figure 17,- Continued. 




0 40 80 120 160 200 

Degrees of freedom 

(f) Convergence for the shear buckling of isotropic clamped 5x1 panels 
\y/\x = = 1.2. From reference 1, Nxy ~ 9.3. 

Figure 17,- Concluded. 







Figure 19, - Variation of i/Ax with trigonometric wavelength parameter 



